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INTRODUCTION 


I  . 


There  are  many  theories  to  aid  the  airfoil  design 
process  by  computational  methods,  because  of  the  desire  to 
reduce  the  number  and  cost  of  wind  tunnel  tests.  A  still 
largely  unresolved  question  is  the  problem  of  flow  sepa¬ 
ration.  Eecause  the  classical  boundary  layer  approximation 
cannot  be  applied  to  separated  flow  calculations,  engineers 
have  tried  to  overcome  this  limitation  by  developing 
v i sccus / inv i sc  id  interaction  approaches  or  to  develop  direct 
solutions  of  the  Navier-5tok.es  equation. 

The  purpose  of  this  thesis  is  to  demonstrate  the 
capability  of  the  v i scous / inv i sc  id  interaction  method  by 
applying  Cebeci's  interactive  computer  program  to  a  single 
airfoil  ( F7  63-137)  at  three  low  Reynolds  numbers  and  by 
comparing  the  results  with  experimental  data. 

Chapter  II  explains  the  boundary  layer  theory.  The 
boundary  layer  equations  are  derived  and  the  turbulence 
models  are  introduced.  Also,  this  chapter  includes  the 
prediction  of  transition  boundary  layer  calculations  and 
flow  separation. 

Chapter  III  introduces  the  interaction  methods. 

Three  weak  interaction  methods  are  explained  briefly  and  the 
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s imul taneous  method  is  presented  as  a  strong  interaction 
method . 

Chapter  IV  describes  Cebeci's  interactive  computer 
program.  Input/Output  data  description  and  JCL  files  are 
included.  Also,  the  results  of  the  application  of  Cebeci's 
program  are  discussed. 


1 


II .  BOUNDARY  LAYER  THEORY 


A.  DERIVATION  OF  BOUNDARY  LAYER  EQUATIONS 

Generally,  the  thickness  of  the  boundary  layer  increases 


with 

viscosity,  or  it 

is  possible 

to  state 

that 

it  decreases 

with 

v  iscos i ty ,  or  it 

is  possible 

to  state 

that 

it  decreases 

as  the  Reynolds  number  increases.  From  exact  solutions  of 
the  Na v ier -Stokes  equations,  it  was  seen  that  the 
boundary-layer  thickness  {(f)  is  proportional  to  the  square 
root  of  kinematic  viscosity  ())); 

cf{x) 

where  x  is  the  distance  from  the  leading  edge  of  a  flat 
plate.  Using  the  local  Reynolds  number 

Re  =  U  x/V, 

For  simplicity,  assume  a  two-dimensional,  steady 
constant  -  property  flow  without  body  forces  and  leave  the 
stresses  unspecified  so  that  the  results  apply  to  laminar  or 
turbulent  flow.  If  simplifications  are  to  be  introduced 
into  Navier -Stokes  equations,  it  must  be  assumed  that  the 
boundary  layer  thickness  is  very  small  compared  with  a 
representative  linear  dimension  (L)  of  the  body,  re.  <f  1  ■'  L. 
In  this  way,  the  solutions  obtained  from  the  boundary  layer 


equations  are  asymptotic  and  apply  to  very  large  Reynolds 
numbers . 

From  the  Navier -Stokes  equations. 


x -direction : 

-  1 

3  P 

+  l 

3^xx 

+  1 

3<Txv 

=  u  3  u  + 

v  3  u 

P 

3  X 

P 

3  x 

/> 

3  y 

3  x 

ay 

y-direction : 

-  1 

3  P 

+  l 

3  6~y  v 

+  1 

3<Txy 

-  u  3  v  + 

v  3  v 

P 

3  y 

P 

3  y 

p 

3  x 

d  x 

ay 

and  the  continuity  equation  is: 

3  u  +  3  v  =0  (2.3) 

3  x  ay 

Inserting  the  "typical"  ( order -o f -magnitude }  values,  replace 
the  dependent  variables  as  follows; 

u  <^>  ,  3  u  ue  3  u  ue 

e  - - —  ,  -  Co  - ~ 

ay  f  ax 

Then  Eq .  (2.1)  can  be  expressed  as 

-  1_  3  P  +  1_  3  <Txx  f  l_  dlTxv  =  u  3  u  +  v  3  u 
/°ax^3x  /°  ay  ax  ay 

u«i*  Q~xx//°  G'xy  //°  ue*  ue^  (2.4) 

11  1  1  1 

where  the  typical  values  were  written  below  the  terms  to 
which  they  correspond. 

Considering  turbulent  flow,  all  stresses  are  of  the 

same  order  ie.  a  general  stress  must  be  of  order 

/°  Ue  / /l,  then  Eq .  (2.1)  becomes 

-  1  d  P  *  1_  3<Txy  1  +  0(  <£_\  ~  u  3  u  *  v  3  u  (2.5) 

P  d  x.  P  3y  '  1  '  J  3  x  ay 


where  0(i/l)  indicates  a  quantity  of  the  order  of  magnitude 

of  in. 

Considering  a  laminar  flow  of  Newtonian  viscous 

fluid, 

(T  xx  =  2 yU  d  u  (T  x y  /  3  u  +  3  y\  (T  y  y  =  2//  3  v 

3  x'  (ay  3x  /  ^y 

In  the  (T xy  term,  dv/ax  is  smaller  than  <Ju/ay. 

Therefore,  Eq .  (2.5)  becomes 

1  c 

-  1_  3  P  4-  y  _d  u  1  +  0  ^ £_  ]  =  u  a  u  +  v  3  u 

/°  3  x  3  /*  |_  ’  1  '  J  3  x  3y 

Similarly,  Eq .  (2.2)  can  be  written  as: 

-  1_  3  P  +  1_  3(Tyy  +  3(Txy  =  u  3  v  +  v  3  v 

/°  a  y  i°  d  y  P  3x  ax  ay 

1  3  P  P  ue  .  /  lie  Ue<T  \  u//  ue  <f  (2.8) 

/°  ay  i/  M  1/  ’  /J  /  i*  I* 

If  we  write  all  the  viscous  terms  together, 

1  3  P  ugV  /  U*  l,(£_)X  u//  u/ /  (2.9) 

P  3  y  1*  ue  1  1  «//  L  ll  /  J  1A  1-* 

It  is  known  that  ( c/7  1  )  <^»  i//ue  1  is  laminar  flow,  so  that  the 

viscous  terms  are  also  of  order  Ug  <f /I  ie.  3P/^y  is  of 

order  UgJ/1  ,  but  the  pressure  difference  between  y=0  and 

y=  i  is  of  order  p  u/ <T /l*and  the  difference  in  dP/dx  will  be 

negligible  compared  to  the  external  stream  dynamic  pressure, 

X 

1/2  />  u e  .  Therefore,  for  practical  purposes, 

3  P  =0  (2.10) 

3 


(2.6) 


(2.7) 


For  this  case,  since  changes  in  P  must  be  of  the  same  order 

as  changes  in  y,  the  pressure  does  not  change  significantly 

through  the  boundary  layer. 

Thus,  the  entire  equation  of  motion  in  the 

y-direction  may  be  dropped  from  further  considerations.  In 

this  way,  the  following  simplified  equations  are  left  for 

the  analysis  of  a  boundary  layer: 

<3  u  +  £ v  =0  (2.3) 

d  x  ay 

-  i  d  p  +  ^  u  -  u  a  u  +  u  a  u  (2.11) 

!°  ax  3  y1  d  x  d  y 

3  P  =  0  (2.10) 

a  y 

These  relations  are  known  as  Prandtl's  boundary  layer 
equations.  Unless  one  encounters  very  high  Mach  numbers, 
the  aocve  orders  of  magnitude  are  not  changed  when 
compressibility  effects  are  considered. 

B.  LAMINAR  AND  TURBULENT  BOUNDARY  LAYER 

The  low  viscosity  fluid  flow  past  solid  bodies  should  be 
considered  as  two  regions.  One  is  the  thin  region  near  the 
boundary  in  which  the  effects  of  viscosity  are  concentrated 
and  the  other  is  the  region  further  away  from  the  wall  in 
which  the  influence  of  viscosity  is  so  small  that  it  can  be 
neglected.  Thus,  it  can  be  stated  that  all  viscous  effects 
are  concentrated  in  the  thin  region  which  is  known  as  the 
boundary  layer.  This  boundary-layer  type  behavior  requires 


15 


high  Reynolds  numbers.  Generally,  the  thickness  of  the 
boundary  layer  (/)  is  defined  as  that  distance  from  the  wall 
where  the  velocity  (u)  differs  by  1%  from  that  (Ue) 
calculated  by  the  ideal  flow  analysis. 

Consider  a  constant-property,  steady,  two-dimensional 
flow  past  a  flat  plate.  If  u/Ue  is  plotted  against  a 
dimensionless  y-coordinate,  T  =  (Ue/vx)  /  ,  the  velocity 
profiles  are  geometrically  similar  and  reduce  to  a  single 
curve  for  a  laminar  boundary  layer  flow.  This  is  well  known 
as  the  Blasius  profile.  The  geometrical  similarity  is 
maintained  regardless  of  the  Reynolds  number  of  the  flow  or 
of  the  local  skin  friction.  For  a  turbulent  boundary  layer 
flow,  since  the  viscous-dependent  part  and  the  remaining 
Reynolds-stress-dependent  part  of  the  profile  require 
different  length  scaling  parameters,  there  is  no  choice  of 
dimensionless  y-coordinate  that  leads  to  the  collapse  of  the 
complete  set  of  velocity  profiles  into  a  single  curve. 

The  conspicuous  difference  in  profile  shape  between 
laminar  and  turbulent  shear  layers  can  be  found  in  the  wall 
flows.  Because  of  the  constraint  on  eddy  size  in  wall 
flows,  the  efficiency  of  turbulent  momentum  transfer 
decreases  rapidly  near  the  wall.  But,  the  efficiency  of 
viscous  momentum  transfer  is  not  dependent  on  y  distance  in 
the  flow  which  has  no  heat  transfer. 
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TURBULENCE  MODELS 


The  unsteady  continuity  and  Nav ier -Stokes  Equations  are 
valid  in  both  laminar  and  turbulent  flows.  Eut,  it  is  too 
difficult  to  deal  with  the  instantaneous  properties  in 
turbulent  flow  because  the  turbulent  flow  has  a  complex 
time-dependent  behavior.  Thus,  the  following  turbulence 
models  are  used  to  make  the  analysis  of  turbulent  flow  more 
convenient . 

1 .  Prandtl 's  M i x i ng -Length 

Consider  two  adjacent  stream  layers  of  fluid  which 
move  with  different  velocities.  If  a  particle  of  fluid 
moves  from  one  layer  to  the  other,  a  momentum  change  occurs 
between  two  layers.  The  fast  particles  which  enter  the  slow 
layer  make  it  faster  and  the  slow  particles  which  enter  the 
faster  layer  impose  a  drag  on  it. 

The  mean  velocity  of  a  stream  layer  is  u,  and  that 
of  the  other  is  u  +  1  d  u/^y  where  1  is  the  distance  between 
two  layers.  Also,  the  fluctuating  velocity  in  the 
x-direction  is  u',  and  that  in  the  y-direction  is  v'. 

Prandtl  assumed  that  the  turbulent  fluctuations  are  due  to 
the  difference  in  the  mean  velocities  of  the  two  layer  3  .  jO 
u '  =  1  a  u/ay  ie.  the  fluctuating  velocity  in  the  x-direction 
is  of  the  order  of  the  difference  in  the  mean  velocities  of 
two  layers  which  have  a  distance  1,  where  1  is  the  mixing 
length.  Prandtl  also  assumed  that  all  components  of 


fluctuating  velocity  at  a  given  point  are  of  the  same  order 
of  magnitude.  Thus,  v'  can  be  defined  as  v*  =  kl  au/ay 
where  k  is  a  constant. 

The  turbulent  shear  stress  due  to  momentum  exchange 
between  two  layers  is  the  rate  of  momentum  transfer  per  unit 
area.  Then  the  mean  turbulent  shear  stress  on  the  fluid  is 
=  -  fi  u  '  v  '  where  u'v'  °*|l  du/«jy|  |  kl  9u/ay|  .  Since  the 
values  of  1  and  k  are  unknown,  combine  these  two  unknowns. 

X,  _  *  _ 

Then,  |Ju/ay|au/ay  where  1  is  called  the  mixing 

length. 


2 .  Cebeci -Smith  Model 

With  the  eddy  viscosity  concept,  the  momentum 
equation  for  2-dimensional  laminar  and  turbulent  boundary 
layers  can  be  written  as: 


(  b  f  "  )  '  +  m  ♦  1  f  f  "  +  m  fl 

.  -  (  f  '  f  ] 

=  x/f'Al' 

-  \ 

2  L 

J 

'  d  x 

d  x  ) 

where  x  is  the  transformed  x-variable 
£  =  ^ m/ V 

Id  -  (  1  +  t )  (  1  +  £  ) 

Let  the  turbulent  boundary  layer  be  a  composite 
layer  consisting  of  inner  and  outer  regions.  Then,  the 
eddy -v iscos  i  ty  formula  for  the  inner  region  is: 


(  >;  =  ll 


(0  i  y  i  yc  ) 


where  1  is  the  mixing  length  for  2 -d imens ional  flow 


is  an  intermi t tency  factor  for  a  flat  plate. 

R*  =  Up  xtr/))  is  the  transition  Reynolds  number. 

The  eddy -v i sco s i ty  formula  for  the  outer  region  is: 

I  co 

(£n)t  -  c/  l(up-u)dy  f  (y  £  y  <_<f ) 

’o 

where  Rg  >_  5000,  the  universal  constant  © <  =  0.0168. 

R$  <  5000,  oi  varies  with  Rg  according  to  the 
the  empirical  formula 

oL  -  0.0168  1  .  55 

1  +  A 

r  “I 

A  =  0.55  [_1  -  exp  (0.243  k  -  0  .  298  K  )  J 

K  =  Re  /  4  25  -  1 

But  this  model  is  not  used  in  Cebeci's  interactive  computer 
program  which  will  be  presented  in  Chapter  IV. 

D.  TRANSITION 

The  boundary  layer  with  a  finite  thickness  starts  out  as 
laminar  at  first  in  the  flow  past  an  airfoil.  However,  the 
boundary  layer  becomes  unstable  and  all  small  disturbances 
begin  transition  to  the  erratically  unsteady  condition  which 
is  known  as  turbulence. 

In  the  boundary  layer  on  blunt  bodies,  transition  makes 
the  point  of  separation  move  downstream  which  decreases  the 
width  of  the  wake.  There  is  an  abrupt  change  in  the  drag 
curve  of  a  sphere.  This  change  is  due  to  a  boundary  layer 
effect  and  is  also  one  of  the  transition  phenomena. 


The  process  of  transition  on  a  flat  plate  at  zero 
incidence  shows  a  sudden  increase  in  the  boundary  layer 
thickness.  Furthermore,  transition  involves  a  great  change 
in  the  shape  of  the  velocity  distribution  and  a  large 
decrease  in  the  ratio  of  the  displacement  thickness  to  the 
momentum  thickness.  Also,  it  causes  a  large  change  in  the 
skin  friction. 

As  described  in  the  above,  the  transition  from  laminar 
to  turbulent  flow  plays  a  very  important  role.  Since  the 
transition  is  not  an  instantaneous  process  and  the  flow  is 
intermittently  laminar  and  turbulent  over  a  certain  length 
of  the  airfoil,  a  certain  point  where  transition  takes  place 
cannot  be  assiyned.  At  present,  there  are  no  exact  methods 
to  calculate  the  flow  within  the  transitional  region. 
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Nowadays,  however,  two  methods  -  Michel's  method  and  the  e 
method  -  are  used  to  predict  the  transition  empirically. 

1  .  M iche  1  '  s  Method 

Michel  investigated  many  kinds  of  data  for 
incompressible  flows  without  heat  transfer  and  found  this 
empirical  correlation, 

-i  0.4-6 

Rj  *  1.174  1  ♦  22400  R*fr 

L  R  x  - 

where  Rg  =  Ue  d/V,  -  Ue  x/v>  . 

Since  the  momentum  thickness  grows  more  rapidly  when 
the  pressure  gradient  is  positive,  Michel's  equation 
involves  the  effect  of  pressure  gradient.  However,  the 


0 


effect  of  surface  roughness,  which  is  also  very  important, 
is  not  included. 

2  •  The  e^*  Method 
<? 

This  e  method  is  of  spatial  amplification  theory 
based  on  the  solution  of  the  Orr -Sommer  feld  equation  which 
forms  the  point  of  departure  for  the  stability  theory  of 
laminar  flows. 


gf  =  a  =  esp  (  (<* i  - 

i  d^r] 

dX 

9  A  XA 

d  fir  ' 

where  X  -  x//*,  ,  4/=^;  <f  * 

g'(x,y,z,t)  is  a  typical  disturbance  function. 

cJ  =  2  7T  (  A  is  the  wavelength  of  the  disturbance)  . 

A 

^r  is  the  circular  frequency  of  the  partial  oscillation 

^r  is  the  amplification  factor  which  determines  the 
degree  of  amplification 

The  basic  assumption  is  that  transition  begins  at 
the  point  where  a  small  disturbance  introduced  at  the 

y 

critical  Reynolds  number  is  amplified  by  a  factor  of  e 

E.  LAMINAR  BOUNDARY  LAYER  CALCULATIONS 
1 .  Similar  Solutions 

a.  Blasius  Solution  For  a  Flat  Plate 

Assume  -  a  flat  plate  at  zero  angle  of  attack 
-  steady,  2-D,  laminar,  incompressible 


-  constant  visccsity 

-  negligible  body  force 


Since  the  pressure  along  the  plate  is  constant, 
there  is  no  pressure  gradient.  The  simplified  Navier  Stoke 
Equa t ion  is: 

JL 

u  d  I]  +  V  a  u  =  y  iJL  (2.12) 

a  x  a  y 

(B,C)  y=0;  u  =  0 ,  v  =  0 

y  -  oo  ;  u  =  U» 

To  transform  from  the  partial  differential 

equation  to  an  ordinary  differential  equation,  define  the 

following  tr a  ns f or ma t ion  parameter: 

7  ~  y  where  7  =  f(x,y) 

J]/x 


4*  ~  Jl?  xL'w  f  where  f  =  f(^)  only 

The  stream  function  was  defined  to  satisfy  the 
continuity  equation. 


ay 


a  x 


a  x 


'xU,  d  £ 
d  7 

if  = 

ay 

U^,  £',  where 

f '  =  df  (2.13) 

d  7 

df  '  i? 
d  "7  ax 

=  u. 

£" y  [¥ 

1 

2x**  _ 

U»f"  ^ 
2x  L 

(2.14 

Jiy  xU* 

,  d  f 

*7  -Jvu» 

(1/2  x 

)  f 

d  7 

ax 

1/2  ~V Ua>  (  f  '  7  -f  ) 


(2.15 


_A i  =  uj 

U*  df" 

=  Ui  f (2.17) 

y*  J 

Wx  d?  ay 

^  x 

Substitute 

Eqs .  (2.13)  ' 

(2.17)  into  Eq .  (2.12),  then 

f  f  »  + 

2  f "  '  =  0 

:  E .  c ;  7 

=  C;  f'  =  0, 

£  -  0 

* 

-co  ;  f'  =  1 

The  s 

olution  of  thi 

s  Blasius  equation  is  presented  in 

" Boundary 

layer  Theory" 

by  Schlichting.  From  the 

transfer ma 

t  i  o  n  relation. 

<f  = 

5 . C  =  5.0 

<f  =  1.7208 

X 

JV o  x / y  J  R  x 

x  /r7 

/  V  )  = 

0 . 8604 

9  =  0.664 

\  U  Iq 

IK 

x  iK 

II 

0.332  /(Ja 

ly« 

\ 

i  l/  X 

^  r  _ 

1  ’ 

0 . 664 

~FT 

V 

Fal kner -Skan 

Method 

The  Falkner-Skan  transformation  is  for  2-D, 
ixisymmetric  laminar  flow.  The  simplified  Mav ier -Stokes 
equation  is; 


u  i  u  +  v  d  u 
ix  d  y 


-  1  P  +  j)  d  u 

(°  *  X  <3  y-1 


'  B,  C )  y  =  0 ;  u  =  0 ,  v  =  0 

y  =  oe  ;  u  =  'J  (  x  ) 


Take  the  same  as  Blasius'  but  different  with  a  function 
£  +  £(x,7)  and  follow  the  same  procedure  as  Blasius'  using 


f  d  x  dx,  then  the  transformed  equation  is: 

V"  +  m  +  1  f  f "  +  m  -  :  f  '  \  J  -  x  ^  f  •  J_t '  f ’’ _±f_  \ 

2  d  x  ^  x  / 


where  m  is  a  dimensionless  pressure  gradient  parameter 
de  f i ned  by : 


x  do 


TJ  dx 

:E.C)  ’{  =  3;  f  =  0,  f  =  0 

7  »;  f  ’  =  1 

For  similarity  in  2-3,  laminar  flow,  assume  f  is  a  function 
of  ^  only.  Then, 
f  +  m  +  1  f  f  ”  +  m 


■n- 


(2.18 


(B.C;  ^  -0;  f*  =  constant,  f 

1  =oo 


ft  -  7 


The  fact  that  m  is  a  constant  leads  tc  : 


c  x 


m 


where  c  is  also  a  constant. 

In  the  case  of  m  -  0,  ie.  U  is  a  constant, 

Eq .  .3.13)  reduces  to  the  Elasius  equation.  If  m  -  1  which 

means  a  2-3  stagnation  flow,  Eq .  (2.18)  becomes  the  riiemens 

e  q  u a  t .on, 

f +  f  f ''  -  (  f  '  )x  +  1  =  0 

Some  solutions  of  the  Faik.ner-Sk.an  equation  for  various 
values  of  m  are  presented  by  Cebeci  and  Bradshaw. 

2 .  Integral  Methods 


ntegral  Momentum  Equation 


a  a  u 
a  x 


(2.19) 


For  steady,  2-D  and  incompressible  flow, 

X 

+  v  a  u  =  -  1  a  p  +  y  a  u 

ay  P  dx  ay1 

a  a  +  a  y  =o  (2.20) 

ax  ay 

From  Eq.  (2.20),  u  /  a  u  +  a  y  \  -  0 

V  a  x  ay/ 

At  the  edge  of  the  boundary  layer,  U  (x)  a  U  ( x )  -  -  l_  a_P 

ax  /’ax 


Then,  Eq .  (2.19)  becomes: 

_a_u A  +  a  ( i j  y ;  =  u  ( x )  d  u  (  x  )  +  a^g 

ax  ay  dx  ay1 

Integrate  this  equation  with  respect  to  y,  from  y 


y  =  /  ,  using 

a  u 

ay 

1  i  X  1 

a  u  dy 

U(  x) 

a  u  dv 

-  |  U(x)  dU ( x  )  d 

o 

«V/ 

X 

1 

0 

a  x 

'o 

dx 

r  * 

(S, 

Also  we  know. 

/  = 

I1  " 

u  ) 

dy 

u 

U(x)  / 

e  - 

> 

(S 

u 

f  1  -  . 

u 

jdy 

0  U(x) 

1 

U  (  x  ) 

Substitute  Eq . 

(2.22)  into 

Eq.  (2 

.  21 ) , 

then. 

dUa(  x  )$ 

+  / 

U  ( x ) 

dU  (  x  ) 

— 

r 

dx 

dx 

0  to 


f  (2.21) 


(2.22) 


This  equation  is  known  as  the  momentum- i ntegr a  1  equation  of 
boundary  layer  theory,  or  as  von  KArmAn's  integral  equation, 
b.  Pohlhausen's  Method 

Assume  a  polynomial  of  the  fourth  degree  for  the 


velocity  function 


where  /  is  the  dimensionless  distance  from  the  wall, 

A  =  y //  and  0  <  A  <  1 

(B.C)  ^  =  0;  f  =  0, 

*  =  1;  f  =  1,  f  =  0,  f"  =  0 

Also  U(x)  d  U  ( x )  +  \)  <?_  u  =  0  when  ^  =  0 

dx  *  yx 

_/_*  d  U  ( x )  =  -  2b  =  A 

y  dx 

where  is  the  shape  factor . 

then,  f  (-*)  =  F  ( A  )  +  A  GU) 

3  4 

where  F(>)  =  2  ^ -  2  A  +  A 
GO  )  =  1/6  A  (1  -b)3 


Thus  the  boundary  layer  parameters  S  ,  6  and  Tw  can  be 
determined,  if  the  velocity  profile  is  known,  as  follows: 


6  10  120 

JL  =  _ L_  /  1Z_  A  -  A  ) 

/  63  \  5  15  144  ' 

Tw  =  ^UJjLL  (  2  +  A  ) 

c.  Thwaites'  Method 

The  integral  momentum  equation  can  be  written 

as : 

d0  +  6  (H  +2)  dU(  x )  =  CJ_  (2.23 

dx  U  ( x )  dx  2 

(B.C)  y  =  0 ;  /u  *  -  U(x)  K  ,  *u  =  U(  x )  L  (  K ) 


where  H  = 

0/  =  Tw  =  ^  L 

2  f>  U*(x)  flUIxl 

K  =  A*  dU ( x ) 

V  dx 

then  Eq .  (2.23)  can  be  reduced  again. 

U( x)  d^L  »  2  |  L(K)  -  [H(K)  +  2]  =  F ( K )  (2.24) 

y  dx  1  1 

Thwaites  writes  an  expression  for  F(K) 

F ( K )  =  0 . 45  -  6K 

Then  we  can  write  Eq .  (2.24)  as 

X  I*  S 

U( x)0  -  0 . 45  U  (x)  dx 

])  UJ(x)  J0 

If  0  is  calculated  for  a  given  ex  ter nal -ve 1 oc i ty 

distribution.  H  and  C  ^  can  be  determined  with  the  following 

r elat ions . 

For  0  i  K  i  0 . 1 

L  =  0.22  +  1 . 57K  -  1 . 8K 

H  =  2.61  -  3 . 75K  -  5 . 24K 

For  -0.1  IK  <.0 

L  =  0.22  +  1.402K  +  0 . 016K 

0.107  +  K 

H  =  0.0731  +  2.088 

0.14  ♦  K 

3 .  Fi nl te -Pi f fer ence  Methods 

The  finite  difference  methods  are  the  most  flexible 
practical  and  efficient  methods  to  solve  the  boundary  layer 


equations . 


A  prime  denotes  differentiation  with  respect  to  7  . 

The  boundary  conditions  are 

7  =  0;  f(|,  0)  =  0,  u(  ^  ,  0)  =0  (2.26) 

7=7  ;  u  (  %  '  la  ]  =  1 

We  denote  the  net  points  shown  in  Figure  2.1  by 

K  =  0;  $  =  *  +  K„ 


?.  - u;  ?/: 

7*> 


+  h; 


n  =  1 ,  2 ,  .  .  .  ,  n 
j  =  1 ,  2,  .  .  .  ,  J 


(2.27) 


Here  n  and  j  indicate  sequence  numbers. 

With  g^  =  g (  £  n,  7  j )  denoting  the  value  of  any 
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quantity  at  the  mesh  point  (  ^ n,  7j),  centered  quantities 


can  be  written  as  two-point  averages: 

p'W  n  n- 1  ~ 

f  =  1  /  2  (  i  +  *  ) ,  1  /  2  <  7j  +  ij-1)  (2.28a) 

an'%  r  fl'i  n  n  n 

<i  *  l/2(  £  ♦  tj).  %tJA  l/2(^  +  ^.()  (2.28b) 

The  finite  difference  approximations  of  Eqs .  (2.25a) 

and  (2.25b)  for  the  mid  point  of  the  segment  P(  ,  are: 

(2.29a) 

( 2 . 29b) 


fl 

f  .r 

n 

f  i-t 

n 

=  u  .  .. 

*-v± 

n 

n 

h> 

n 

U  J, 

u  j  -• 

1 

ii 

h> 

Eq .  (2.25c)  can  be  approximated  similarly  by 

centering  about  the  mid  point  of  the  rectangle  P1,P2,P3,P4 

1.  Centering  Eq .  (2.25c)  about  the  point 

.n-'A 

1  %  ■  n  > 

Let  the  left-hand  side  of  Eq.  (2.25c)  be  L  and  use 
Eq.  (2.28b) 


29 


b. 


Newton's  Method 


For  simplicity,  let  (  f  ?  ,  u-  ,  v^  )  be 


£j ,  uy ,  v-  )  at  i  = $  "  . 


Then 

Eq. 

(2. 

29)  and  2 

.  31  ) 

can  be 

written  as 

- 

c  •  _ 

‘■i- » 

h} 

2 

( ua-  + 

u*-f  )  =  0 

(  2  . 

- 

U3-l  “ 

M 

2 

( va  + 

v/- ,  5  =  0 

(2. 

.  1 

h  - 

3 

1 

(  b 

3  V  -  b3- 

-1  v 

/- 1  )  + 

<k,  (  fv  ) 

(  u  )  y 

4* 

n  - ) 

*(V*  f 

9-* 

n  - 1 

~  £t-Jl 

n  -1 

5  =  R?-*. 

(  2  . 

Here 

the 

unk 

now ns  are 

on 

the  left 

-hand  side  and 

n  -  1 

R  3  ' 

involves  only  known  quantities. 


Now,  Newton's  Method  is  applied  to  turn 


Eq .  ( 2 . 3  3  )  into  a 

linear 

system 

complemented 

by  boundary 

conditions  Eq . ( 2 . 

32)  and 

initial  values 

rt> 

fl 

0 

(a 

= 

0 

(*) 

Vo  = 

Vo"" 

(»■>  *  *  * 

f  •  -  i  . 

3  3 

(*) 

*9  = 

n'1 

U  • 

3 

ii 

5,  'r\> 
> 

n  -» 

V 

(1  <  j  '  J 

<•> 

J  or 

W 

UJ  = 

1 

(J) 

n  -* 

vj 

The  superscripts  in  parenthesis  represent  the  iteration 


number  as  follows: 


to 

' 


i  =  0, 


2, 


f*5  (■  (141)  r  (l>  (l*<)  (l>  (■  «i> 

t  +  /£,-,  =  u#-  +  /u^-  ,  V.  =  vy  ♦  /v#. 

where  /f  <<  f,  /u  <<  u ,  /  v  <<  v. 


(i*<)  h> 

f,  £ 


Replace  fi,  Uj,  v^  in  Eq .  (2.33)  with  these  expressions 

Then  Eq .  ( 2 . 3  3a )  becomes : 


/ft  -  If-  - 

.  ,  Sty  .  "i 

u  T 

L;u/  ^u^-'  J 

=  ( r,  )y  (2.34a) 

where  (  r  ,  )  -  =  f  y_ , 
From  Eq  .  (2.33b), 

til  ,  ll>  to 


tt  +  h3  u#'* 


From  Eq .  (2.33c) 

h  ^  (  bj  /  vy 

*  -  • 

+  (  V V 


.  -  /  V,  ~ 

hi 

2 

in  r  m  tn  -  i/i 

L  )  *  '  u)  *  v.  *  *  V. 

¥  [/u>  * 

]  =  l- 

ID 

-  U  .■ 

3->  )-< 

/ 1) 

-  Uy 

+  h3  V*-Ji 

to  r  II) 

by-,  /  v,., 

)  ♦ 

,  f  " '  r  a  "  > 

<*,  0  (  fv  )y.  ^  -  c*,/  (  u  )y,  yx 

r  *'» 

/£/-x  ■ 

rt  - » 

r  Ml 

'  VJ-  4  J 

2  .  3  Ab  ) 


n  - 1 

--  R,*.* 


(n  Cl 

h  y  (by  Vy  -  by.t  vy.,  '  *  U,  (  f  r 


•  -  <  >;:*  iA  -  K-l  <  ) ' 

her.  [fjVv?  •v;-/fy.  f" 

r  II I  -  in 

'  (  J  V-y*  -  U,  '  u#  +  V<  >  U  y.. 


,  J  <1  » 

<V  »  '  u  *  , 

<r/ 


/*  t*  > 

*  f  J  *  '4 

r  (i) 

^  v  • 


1/2  (  /  fy  ♦  /  f  y_  ,  ) 

1/2  (  ^  v  '■  ’  ♦  /  v  ) 

?  V- ' 


then  (s,)y/vy  +  (saLVv.„  +  ( 5j  ).  ^  f  +  (s„'-/f3 
♦  <s<  )y  /uy.,  =  (ra)y  *  * 


*  c\r  j  ■ 

(2  ?  A . 


(  r  i )  y 

*  -  1 

=  Rj-4  “ 

- 

* 

<=,  v 

"  1  IM 

=  hy  by  ♦ 

<Si’/ 

- 1  ( 1  / 

=  "  hi  *V- 

<S,V 

,  <•> 
=  Vy 

♦  Ot ,  ,  f  V 


*  -t  _  fn 

V* 

_  M  > 


.  n  ’  •  I  ; 

9-  *  V-  >i 


*  -  I 

V  - 

3  -  '4 


3  2 


(s4)3-  = 

e^i 

2 

(I  r 

V.  + 

-  v' 
2  3  ' 

II 

** 

l/) 

- 

,  »•' 

ts(>-  = 

- 

o/A  u  |  ’  , 

the  boundary  conditions  are: 

/ 10  =  0  /u0  =  0  f  Uj  =  0  (2.35) 


| 


c.  Block  Elimination  Method 

Eqs.(2.34)  and  (2.35)  are  the  linearized 
difference  equations  of  the  momentum  equation  for  external 
flows  which  has  a  block  tridiagonal  structure.  We  can  solve 
these  equations  by  means  of  the  Block  Elimination  Method  as 
discussed  by  Keller  (1974). 

Let  us  define  the  three-dimensional  vectors,  / - 

d 


and  r  •,  to  express  the  system  in  ma t r  i  x -vec t or  form, 
<f 


where 


A q  Go 

B#  A  f  C | 


V 

t 

'  r»  " 

b3  A/  CJ  _ 

/  = 

4 

f?  = 

r, 

0 

- 

/j- 

^ J'’ 

Bj-i  Aj.,  Cj.t 

B  7  A  j 

— 

/  7 

_rJ, 

Let  as  factorize  the  matrix  ft  to  solve  Eq .  (2.38) 

(S  =  xy 


where 

[i 


P> 


X  ,  I 


1 


*3-i  I 


XJ  1 


y.  c* 


(2. 39) 


V  '  c, - 

i 


h., CJ- 


Here  I  is  the  3x3  identity  matrix. 

x„-  and  y„-  are  also  3x3  matrices. 
According  to  Eq.(2.39),  we  can  find 
Yo  ~  A  0 

Xj  yy =  B^  (  j  =  1,2,...,  J ) 

y}-  =  Ay  -  Xj-  cy.t  (j  =  1,2...,  J) 


(2.40a) 
( 2. 40b) 
( 2 . 40c) 


and  the  matrix  xj  has  the  same  structure  as  that  of  the 
matrix  Bj.  Therefore,  x^  has  the  elements  like  this. 


(  x  „  )  j 
(  x  *,  ) 

0 


i 


(  x  j  j  )  ^  (  x  ,  ^  )  y 

(  X  ^  )  •  (Xj,)/ 

0  0 
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and  y *  can  be  denoted  as: 

4 


^  y  ii  V  ;Yij-  (Yu  y 

(y*.  y  <y«  (y^j  y 

o  -r  -h1+l  /  2 


From  Eq . 

( 2 . 40a  ) 

<  V  ii  '0 

_  t 

—  X 

'i’n  '0  =  C  * 

y ,3  ’a  = 

0 

:y»i  ’0 

_ 

—  <J 

:y~  0  =  1 

y*j  >o  = 

0 

V 

From  Eq . 

2  .  40b) 

(X  1,  ) , 

—  —  1 

(  x  ,x  ) ,  =  -1/2 

h  ,  (  X  ,3  ) 

,  =  0 

‘  X  x,  > , 

=  <s*). 

1  *■  x  u.  ^  1  -  v  ^  1 

+  (x  o>, 

(  xAJ  ) 

From  Eq . 

( 2 . 40c ) 

II 

1  (y 

U  '2  ~  3  _  x/3 

*  2 

V  (y,3  } 

II 

rol?: 

(y^  v = 

(s?  V  i  y. 

u  y  -  (srV"  (xj 

3  V  «y^ 

y=  (s 

Then  we  can  compute  the  elements  of  x  with  Eqs .  (2.3' 


2. 40b; 

i  for  1  i  j  i 

T 

^  / 

x  11  y 

=-r  i  ;Vj’ 

)  - 

3-' 

*  Ja  fii> 

<y*i  ’3'., 

-  <y« 

V->J  ) 

5r  0  1 

L  3 

x  iu  y 

■ "  i  f(^) 

a 

* 

:  x "  v  [ ; 

y.,  - 

h_ 

y  i  3  J-  < 

y,  i\3/ 

Xl3'/ 

=  J2i  t  x  .« 

2 

V  ; 

y/J  V-,  + 

;  x  m  y 

'i’iJ  J-i 

1 

J 

X  11  3 

U) 

r  -« 

II 

<y*i 

J  *  iS«'j 

'7  i3  V- 1 

- 

y*  1 

9 

' s*  y 

iy^  y.r 

(s*  y 

X  »*  V 

=  -4- 1  ^(s 

*  V 

-  (s*y  -  ^ 

X  A,  'y  ( 

y  t)  V- 1 

yi  y. 

c*jV  =  t*»i  V  (^u  V-  +  ix«  V  ;yu  V-« 
y„  -  v y , i  ■  j-i  ^ y ».i  1  ~  '  *  *  3  3- >  n  3-1 


J-l 


Let  y  /  =  z,  then  xz  =  R 
Thus,  z„  =  r0 


(2.41'/ 

(2.42) 


“3  -  L° 

where 


x  •  z  • 

9  3" 


(1  <  j  <  J  ) 


(Z,  )j 

(  za  )y 
(  z  ,  ) !• 


0  i  j  £  J) 


From  Eq .  (2.42)  for  j  =  0 

(z,  >0  =  C  r,  )0  (zj0  -  (£  )„  !zj).  ^  (r,  ), 

and  for  1  <  j  <  J 


(xM 

V ( 2  ■  V- 

-  (  x  (i 

V  (Z/V- 

-  (  X  )/(  Z  J 

v- 

(xi( 

y ( z*  v-< 

-  ( x^ 

V  ^  z-»  V-( 

-  (XAj  )j  (Zj 

V- 

(Z3])  =  (  r3  )j 

From  Eq .  (2.41), 

Vj  /3  ‘  z7  (2.43a) 

yj/,'  =  V  -  =</,•,,  "HU-H  (2.43b) 

Then  the  vectors  can  be  calculated  with  Eq .  (2.43).  The 

three  components  of  ,  for  j  =  0 , 1 , 2 , . . . J - 1 ,  are: 

/  (Jr  =  -  _h|V.  /vv  -  ef 
2 

v,-  -  i  {(y„v[(^v  +  e,(yji)-]  -  <yA,v(z,),'  -  e,(y>,V(y^V  ) 

y* 

[<*.V  -  )■  -  <v^v] 


_J _ 

i  y<« 


where  e,  =  (Zjb  -  o  u^(  +  _h£t.  f  v; 


y^=  |in  v(y'^;j  +  (v~V(y'<  }j 

-  (yit  >/  (y/j  V  +  ;yi3  >/  (y ■.  V 

and  for  j  =  J 

a  j  =  >  2  j  >  j 

6  Vj  =  ej  (y,,  )j  -ej(y„)j _ 

^ y 1 3  'j'Y u  ' j  -  ' y x>  'j 1  y  ii  'j 

'  *j  =  ea  -  (y,j  ).T  /  vy 

(y-  1 J 


whet  e  e  x  = 

-  '  y  u 

;  j/u 

e  3  = 

;  zx  5  J 

-  fy^i 

)jyu 

These  calculations  are  s 

topped 

when 

1  6  v„ 

1  < 

£ 

where  v(o!  is  the  wail  shear  parameter. 
£  is  a  prescribed  value. 


F.  TURBULENT  BOUNDARY  LAYER  CALCULATIONS 

Turbulent  fluid  motion  is  an  irregular  condition  of 
in  which  the  various  quantities  show  a  random  variation 
time  and  space.  Therefore,  turbulence  is  char acter  i  zed 
random  and  chaotic  motion  of  fluid  particles. 

The  velocity  varies  randomly  at  any  point  in  a  turb 
fluid.  The  velocity  components  in  a  three-dimensional 


f  1  o 
wit 

by 


u  1  en 


turbulent  fluid  are: 


where  a,  v, 

u  = 


w 


represent  the  mean  velocities 
1  1  u  dt,  and  so  on. 

tx  -  t,  J. 


u',  v',  w*  represent  the  amount  of  fluctuation  of  the 
instantaneous  velocities  from  the  mean  velocities. 

Obviously,  we  can  see  that 

ftz  — 

1  u1  dt  =  u'  =  o,  etc. 

tx  -  t,  Jit 

Similarly,  the  following  quantities  can  be  defined, 

— 

1  \  u ' *  dt  -  u*  ,  etc. 

tA  -  t,  4, 

,U  - 

1  u'v1  dt  =  u'v’,  etc. 

tA  -  t,  K, 

Even  though  the  flow  is  turbulent,  the  time-dependent 
momentum  equations  are  valid.  However,  it  is  impossible  to 
solve  the  momentum  equations  for  turbulent  flow  because  the 
fluctuations  are  random  and  chaotic.  Thus,  Reynolds 

modified  the  momentum  equations  by  introducing  the  mean 

U. 

value  cind  fluctating  values  of  the  flow  quantities  and  by 
assuming  the  fluctuations  to  be  continuous  functions  of  time 
and  space.  From  the  momentum  equations. 


U  +  U  d  U  +  V  <3  u  +  w  &  u  =  ^_1_  £  P  +  p  (  &  u  +_£_u  +  _£ja  \ 

at  ix  ay  d-  z  f>  d  X  (^X-1  ay1  gzV 


By  multiplying  the  continuity  equations  with  u 


Sum  these  two  equations  and  rearrange,  then 


f^au  - 

u*\  +  3 

1  Pau 

-  uv)  4  3  j 

Pau 

-  uw  ) 

)  <3y 

ay 

1  3z( 

a  z 

/ 

Here  take  the  mean  value  of  each  term  using  u'  =  0, 
3(  u  +  u '  )  =Ju.  +  3u  '  =  a  u 

at  at  3  t  <3  t 

-  _L  *<£  +  P' )  =  zl 

/°  ax  ^  ix 


n  u  j-  u  * ) 


a  ( u  +  u ' )  =  9  a  u 


(u  +  u')X  =  P3u  -  ( u'*’  +  2  uu '  +  u) 

ax 


=  >LiiL 

_  x 


u  -  u 


^  3(  u  +  u  * ) 
a  z 


^  3  u  -  uw  -  u ' w ' 
a  z 


(2.45) 


Substitute  Eq .  (2.45)  into  Eq .  (2.44)  and  rearrange.  Then, 

a  (  ilL  +  u  iiL  +  ^  iA  +  **  ^  \ 

(  \ at  ax  ay  az/ 


/O  (  J-iL  +  u  +  v  +  **  ^  \ 

(  \ at  ax  ay  az/ 

=  -  3 P  ^  P (  3u'V  a u ’ v  1  +  au'w1  j 

ax  ax  ay  a  z  > 


(  2 . 46a  ) 


By  similar  procedures  for  the  other  directions, 
a  (  ay  +  u  3v  +  v  3v  +  w  a  y  N 

'  at  ax  az  / 

=  -  a  p  7  -  a!  a v 1  u 1  +  a  v'a  +  a  v 1  w 1  'i 

ay  '  |  x  ay  d  z  ) 

(  a  w  +  q  a  w  +  v  3  w  +  w  3  w  \ 

•at  ax  ay  3z/ 


( 2 . 46b) 


=  -  jLe 

a  z 


Vf'T"  -fit  d  w  1  u  1  +  i.  w.'v  1  +  -i-W  */  ) 
M  3  X  ay  a  z  / 


( 2 . 46c) 


and  from  the  continuity  equation, 
a  u  +  3v  +  3  w  =  0 


(2.47) 


GW 


In  laminar  flow,  we  have  4  unknowns  and  4  equations. 


But  in  turbulent  flow,  we  have  10  unknowns  and  only  4 
equations  as  shown  at  Eqs.(2.46),  (2.47).  This  is  the 

reason  why  we  need  the  turbulence  models  which  were 
discussed  in  Section  C. 


3.  SEPARATION 

In  some  cases,  the  boundary  layer  thickness  increases 
considerably  in  the  downstream  direction  and  the  flow  in  the 
boundary  layer  reverses  its  direction.  The  change  in 
direction  causes  the  decelerated  fluid  particles  to  move 
outwards,  which  means  that  the  boundary  layer  is  separated 
from  the  wall.  This  phenomenon,  boundary  layer  separation, 
is  always  related  to  the  formation  of  vortices  and  to  large 
energy  losses  in  the  wake  of  the  body. 

Let  us  consider  the  simplified  boundary  layer  equation 
in  order  to  investigate  two-dimensional  separation; 


u  j  u  +  v  £  u 
<^x  ay 


1  d  p  +  u 

f>  **  ayA 


and  define  the  separation  as 


(  £  u 

l  *  y 


Since  u  =  v  =0  at  /  =0, 


The  velocity  profile  in  the  boundary  layer  always  has  an 
inflection  point  in  the  region  of  decelerated  flow  and  the 
velocity  profile  at  separation  point  must  have  an  inflection 
point.  Thus,  separation  can  occur  only  when  the  flow  is 
decelerated.  In  two-dimensional  separation,  a  bubble  of 
fluid  which  has  low  velocity  is  always  formed.  This  bubble 
is  often  unsteady  and  distinguished  by  interior  streamlines 
which  are  closed  loops  or  which  extend  to  infinity. 

There  are  two  severe  problems  in  boundary  layer 
calculations  when  separation  occurs. 

1.  The  Goldstein  singularity  at  the  separation  point  in 
direct  boundary  layer  calculations. 

2.  Numerical  problems  downstream  of  the  separation  point. 

If  a  boundary  condition  prescribes  the  pressure 

gradient,  then  the  boundary  layer  methods  suffer  from  a 
singularity  due  to  separation.  Moreover,  the  singularity 
causes  the  numerical  breakdown  in  direct  boundary  layer 
methods  near  the  separation  point.  This  Goldstein 
singularity  can  be  overcome  by  using  the  displacement 
thickness  or  the  wall  shear  stress,  instead  of  the  pressure 
gradient,  for  the  boundary  condition,  which  makes  it 
possible  to  integrate  the  boundary  layer  equations  through 
the  separation  point.  Also,  the  full  Nav  ier -Stokes 
equations  exhibit  no  singular  behavior. 
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On  the  other  hand,  the  boundary  layer  equations  lead  to 
a  numerical  instability  in  the  region  of  reversed  flow.  The 
FLARE  approximation  is  the  most  common  method  to  overcome 
this  instability.  The  momentum  transport  term 
u  ^u/<?x  is  deleted  where  u  is  smaller  than  zero.  But  the 
accuracy  of  this  approximation  decreases  as  the  region  of 
reversed  flow  increases.  Therefore,  it  is  necessary  to 
introduce  an  upstream  influence  in  that  region.  This  is 
incorporated  by  the  downstream  upstream  iteration  procedure 
(DCJIT)  which  consists  of  a  sequence  of  alternating  up-  and 
downstream  sweeps  with  the  momentum  transport  term  u  au/ax. 


Ill .  INTERACTION  METHOD 


A.  INTRODUCTION 

If  the  flow  remains  attached  and  the  Reynolds  number  is 
high,  the  pressure  distribution  and  the  overall  lift  force 
can  be  obtained  from  a  potential  flow  solution. 

Conventional  boundary  layer  methods  provide  additional 
information  about  the  skin  friction  distribution  and  the 
overall  drag  force.  However,  if  the  flow  separates,  no 
information  is  available  for  regions  downstream  of  the 
separation  point.  Therefore,  the  overall  forces  cannot  be 
obtained . 

For  this  reason,  interaction  methods  are  introduced  to 
overcome  this  problem.  They  provide  a  special  coupling 
between  the  inner  viscous  and  the  outer  inviscid  flows,  and 
can  be  classified  as  follows: 

1.  Weak  interaction  methods 

a.  Direct  method 

b.  Inverse  method 

c.  Semi-inverse  method 

2.  Strong  interaction  method 

The  weak  interactions  provide  only  a  loose  coupling 
between  viscous  and  inviscid  regions,  i.e.  two  different 
regions  are  treated  alternately.  The  viscous  flow  solver 
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deals  with  the  flow  in  the  viscous  region  and  yields  the 
boundary  condition  of  the  inviscid  region,  while  the 
inviscid  flow  solver  deals  with  the  flow  in  the  inviscid 
region  and  yields  the  boundary  condition  of  the  viscous 
region.  The  exchange  of  information  through  the  boundary 
condition  is  slow,  but  regions  such  as  those  near  separation 
and  trailing  edges  require  fast  and  direct  coupling  between 
viscous  and  inviscid  region.  While  the  strong  interaction 
method  treats  displacement  thickness  and  external  velocity 
simultaneously,  the  weak  interaction  methods  process  one  of 
these  quantities  as  input  and  the  other  as  output. 

1.  Inviscid  flow  methods 

a.  Direct  boundary  conditions 

(1)  Prescription  cr  the  airfoil  shape 

(2)  Zero  normal  velocity  at  the  surface 

b.  Inverse  boundary  conditions 
Prescription  of  a  velocity  distribution  for 
the  unknown  airfoil  shape 

2.  Viscous  flow  methods 

a.  Direct  boundary  conditions 

(1)  No  slip  condition  requiring  zero  normal 
and  zero  tangential  velocity  at  the 


sur  face . 


(2)  Prescription  of  the  external  velocity 
i.e.  u-component  of  velocity  at  the 
edge  of  boundary  layer 

b.  Inverse  boundary  conditions 

(1)  No  slip  condition 

(2)  Prescription  of  the  displacement 
thickness 

c.  Boundary  conditions  for  simultaneous 
interact  ion 

(1)  No  slip  condition 

(2)  Prescription  of  a  linear  combination  of 
displacement  thickness  and  external 
velocity 

B  WEAK  INTERACTION  METHODS 
1 .  Direct  Method 

A  direct  inviscid  flow  solver  is  combined  with  a 
direct  viscous  flow  solver  (see  Figure  3.2).  The  boundary 
conditions  are: 

u  (  x ,  o  )  =0  v(x,o)  =  0  (3.2) 

u  (  x ,  y  )  =  Uc(x) 

This  method  is  terminated  at  the  point  of  vanishing 
skin  friction.  In  the  direct  scheme,  the  pressure  is 
calculated  from  the  inviscid  region,  but  the  displacement 
thickness  is  determined  from  the  viscous  region.  Because  o 


this  phenomenon  and  Goldstein  singularity,  this  method  is 
not  appropriate  for  flows  w*th  strong  interference  effects 
between  viscous  and  inviscid  regions. 

2 .  Inverse  Method 

This  method  consists  of  an  inverse  inviscid  and  an 
inverse  viscous  flow  solver  ;  see  Figure  3.31.  The  boundary 
conditions  of  the  inverse  boundary  method  are: 

u  ( x , o  )  =0,  V  (  x  ,  o  )  =0  (3.2; 

4  x ,  y,  )  -  usx,  y e)  [  y«  -/*(x)J 

where  ^  is  the  stream  function. 

The  roles  of  displacement  thickness  and  external 
velocity  distribution  are  exchanged  in  this  method.  The 
troubles  related  to  the  Goldstein  singularity  can  be 
overcome,  but  the  whole  procedure  takes  a  long  time  due  to 
very  slow  convergence.  Thus,  this  method  can  be  applied  tw 
the  regions  of  separated  flow  only  and  needs  severe 
under -relaxation. 

3 .  Semi-inverse  Method 

A  direct  inviscid  flow  solver  is  combined  with  an 
inverse  viscous  flow  solver  (see  Figure  3.4).  The  input  is 
the  displacement  thickness  and  the  output  is  the  external 
velocity  distribution  for  both  solvers.  The  two  external 
velocity  distributions  are  combined  and  then  an  updated 
displacement  thickness  is  obtained  through  a  relaxation 
procedure.  A  formula  for  satisfactory  convergence  is: 
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r  * 

p 

- 

(  X  ) 

=  0  o  1  d  (  X  ) 

1  +  w  / 

h'  evW  -  1  \ 

L  l 

Velffi)  J  _ 

where  w  is  the  relaxation  parameter. 

The  numerical  weaknesses  can  be  overcome,  but  the 
coupling  is  still  loose. 

C .  STRONG  INTERACTION  METHOD 

The  simultaneous  method  solves  the  boundary  layer 
equations  subject  to  an  interaction  law  (see  Figure  3.5). 
Viscous  displacement  effects  are  allowed  to  cause 
substantial  changes  in  the  external  velocity  distribution. 
Since  both  displacement  thickness  and  external  velocity  are 
treated  as  unknowns,  one  more  additional  relation,  the  sc 
called  interaction  law,  is  needed.  Thus  the  boundary 
ns  are: 

*  *  X  ,  O  i  =  v  V  v  X ,  O  )  — 

^  ;x,y>  ;  =  u;x,y€ ;  [  y*  -  o'*  *']  :  - 


~  ^  ■■  *-  = 


X  , 


C*  (  x  '/  *■ 


1  /'  TT  (  — d_  u  ;  i  ,  ye  :  a  ;  i  :  1 _ 

I  d)  L  -*  V 


ci 


'he 


st  equation  represents  the  interaction  law. 

Here,  the  external  velocity  for  the  boundary  conditions 
can  be  written  as: 


U  (  x ,  ye  )  =  u  e  ;  x  ) 


u  o  ;  x ) 


+  & 


b  i  A  / 


(  3  . 


where  Ue  (x)  is  due  to  inviscid  flow  past  the  airfoil 


i x )  is  the  perturbation  velocity  due  to  the 


displacement  effect  of  a  boundary  layer 


To  obtain  41%  ( x ) ,  the  blowing  velocity  concept  is  used.  Le t 
us  consider  an  airfoil  along  which  sources  are  distributed 
(see  Figure  3.1).  This  surface  distribution  of  sources 
represents  the  effect  of  boundary  layers.  There  are  two 
ways  to  predict  the  displacement  effect  of  a  boundary  layer 
on  the  outer  inviscid  flow: 

1.  Compute  the  inviscid  flow  past  a  displacement 
body 

2.  Replace  the  condition  of  zero  normal  velocity  on  the 
surface  with  a  condition  of  prescribed  blowing 
velocity  at  the  surface 


t* 

v  (x,  /  ) 


Figure  3.1  Blowing  Velocity  Concept 

The  streamlines  are  displaced  away  from  the  surface  by 
the  distributed  sources  which  eject  the  fluid  at  the 
surface.  Then  the  virtual  displacement  body  becomes  a 
streamline  and  the  flow  tangency  condition  is: 

v  (  x,  /  )  =  J  </  (3.6) 

Utf  (  x  )  J  x 

From  the  thin  airfoil  approximation,  the  displacement 
thickness  is  assumed  to  be  so  small  that  u-component  of 
velocity  do  not  change  across  the  layer  and  the  airfoil  in 


49 


inverse  Boundary  Ccndit 
Input:  Displacement  thi 

distribution 

Output:  External  velocit 
distribution 


- $ - 

INVISCID  FLOW  SOLVER 


.'Inverse  Boundary  Condi  ti 
Input:  External  velocity 

distribution 

Output:  Displacement  body 
s  tape 


,  * 

c  :  x : 


I NVI SC I D  FLOW  SOLVER 


VISCOUS  FLOW  SOLVER 


(Direct  Boundary 
Conditions) 

Input:  Updated  shape  of 
displacement  body 
Output : External  velocity 
distribution 


(Inverse  Boundary 
Condition) 

Input:  Displacement  thick¬ 
ness  distribution 
Output:  External  velocity 
distribution 


UeV  (*) 


^  CONVERGENCE  ? 
Comparison  of  external 
velocity  distributions 
\Uei  -  UeV  =  0>- 


UelW  \m  U  eV(*) 


RELAXATION 


Update  displacement  thickness 
distribution  thickness  distribution 
and  the  two  external  velocity 
distr  ibutions . 


Figure  3.4  Semi-Inverse  Method 


INVISCID  FLOW  SOLVER 


Air  foil 
Geometry 


(Direct  Boundary  Conditions) 
Input:  Airfoil  shape 
Output:  External  velocity 
Distr ibution 


VISCOUS  FLOW  SOLVER 

_ and  INTERACTION  LAW _ 

(Displacement  thickness  and  external 
velocity  are  treated  as  unknowns) 


Input : 


Output 


Data  obtained  from  previous 
sweep  and  initial  inviscid 
solution 

Displacement  thickness 
distribution  and  external 
velocity  distribution 


i  0  (  X  ) 

RELAXATION  OF 

S*  (  x  )  1 

i 

1 

^^^C^VERGENCE?^-^^ 

1 

i 

S'  Comparison  of  old  and 

No  1 

new  displacement 

S' 

thickness 

Figure  3.5  Simultaneous 

Me  thod 

IV.  DESCRIPTION  OF  CEBECI 1 S  INTERACTIVE  PROGRAM 


A.  INPUT  DESCRIPTION 
1 .  Introduction 

The  Reynolds  number  is  the  most  powerful  parameter 
which  can  affect  the  flow.  Wholly  different  flows  can 
result  from  different  Reynolds  numbers  (e.g.  6.0E  +  06 
against  0.28E  +  06).  Even  though  a  flow  with  high  Reynolds 
number  behaves  well  and  does  not  separate,  it  may  have  an 
extensive  flow  separation  at  a  low  Reynolds  number. 

Numerical  breakdown  in  boundary  layer  computations  is  partly 
due  to  this  extensive  flow  separation  at  low  Reynolds 
numbers.  To  avoid  unrealistically  large  regions  of  flow 
separation,  it  is  necessary  to  make  some  changes  in  the 
turbulence  model  for  transitional  flow.  Reduction  in 
transitional  region  to  a  proper  level  is  a  method  to 
decrease  numerical  problems  in  computation.  By  doing  this, 
more  reasonable  results  can  be  obtained  in  low  Reynolds 
number  flows.  Also,  the  computation  at  low  Reynolds  numbers 
needs  more  iterations  to  converge  than  at  high  Reynolds 
number,  and  the  lift  coefficient  increases  with  Reynolds 
number  at  the  same  angle  of  attack  because  low  Reynolds 
number  flows  exhibit  a  stronger  displacement  effect  than 
high  Reynolds  number  flows. 


55 


The  results  may  be  affected  by  transition  location. 
The  transition  location  can  be  either  fixed  by  the  user  or 
computed  from  the  following  empirical  formula  given  in 
Cebeci  and  Bradshaw  (1977). 

o,  4  6 

R0  =  1.174  ^  1  +  22400  \  Rx  (4.1) 

Rx  / 

Laminar  flow  separation,  however,  may  occur  upstream 

of  the  transition  location  predicted  by  the  above  equation. 

In  this  case,  it  is  assumed  that  the  onset  of  transition 

corresponds  to  the  point  of  laminar  separation,  and  the 

following  message  will  be  issued  in  the  output: 

TRANSITION  LOCATION  HAS  BEEN  RECOMPUTED  AT  THE 
POINT  OF  LAMINAR  SEPARATION 

The  results  of  this  program  agree  well  with  the 

experimental  results  up  to  stall  at  high  Reynolds  number 

flows.  However,  low  Reynolds  number  flows  may  not  agree 

well  or  experience  numerical  breakdown  close  to  stall  with 

the  following  messages: 

at  the  very  top  of  output, 

+  IFY  002  I  STOP  1 

or  at  the  very  bottom  of  output, 

IFY  207  I  VFNTH :  PROGRAM  INTERRUPT  (2)- 
FLO ATI NG- POINT  EXCEPTION  OVERFLOW 

IFY  259  I  STNCT :  /ARG/  =  /argument/, 

(HEX  =  hexadecimal),  APPROACHES  SINGULARITY 
or  MAXIMUM  NUMBER  OF  ITERATIONS  EXCEEDED 

Thus,  the  range  of  angle  of  attack  should  be  considered 


c a refull v. 


If  the  user  requires  many  sweeps  with  IPRNT  =  0  or 
2,  then  the  user  needs  to  choose  class  "G"  or  "J"  for  enough 
running  time  and  add  additional  command  ",  LINES  =  (m)"  just 
after  nnnnP  on  the  following  line  to  make  enough  space  for 
print: 

11*  MAIN  ORG  =  NPGVM1  nnnnP  ,  LINES  =  (m) 

where  nnnn  is  user's  IS  number. 

m  is  the  number  of  output  lines  in  thousand 
(Generally  10 js  enough) 

If  the  user  does  not  do  that,  the  following  error 
message  will  be  issued  at  the  very  top  of  the  output  and 
printing  will  be  stopped  abnormally: 

IAT  1600  JOB  number  (user's  job  number)  LINES  EXCEEDED 
IEF  4501  user's  job  number  -  ABEND  S722  U0000 

2 .  Detailed  Input  Data  Description 

The  input  to  the  computer  program  consists  of  1 

title  line,  3  control  lines  and  airfoil  coordinate.  User 

must  follow  the  data  format  for  specified  column  and  type. 

I  Title  Line  I  FORMAT  (18  A4) 

This  line  provides  any  description  as  desired  with 
any  acceptable  machine  characters. 

I  Control  Line  1  I  FORMAT  (5  15) 


'-fill 


»  *  .  1  »  '  »  *  .  *  >  *  *  * 
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Column 

5 

10 


15 


Name 


Explanat ion 


ITF(l)  Transition  flag  for  the 

lower  surface 

I TF ( 2 )  Transition  flag  for  the 

upper  surface 
=  1  Point  of  transition 
has  to  be  specified 
by  user.  This  option 
shall  be  used  if  ex¬ 
perimental  results 
are  available  or  to 
avoid  oscillations  of 
the  computed  transition 
points.  In  the  latter 
case,  fix  transition  at 
the  most  upstream 
occurrence  of  computed 
tr ans i t ion . 

=  4  Point  of  transition 
will  be  calculated 
according  to  Eq.  C4.1). 
No  input  is  necessary 
for  XTRL  and  XTRU.  The 
computed  point  of 
transition  will  be 
redefined,  if  laminar 
separation  takes  place 
upstream  of  the 
location  predicted  by 
Eq.  (4.1) 

IRSTRT  Boundary  layer  restart 

f  lag . 


Read  Starting  Solution  Store  Final 

_ _ solution  or  Unit  2 


-  0  No  Yes 

=1  Yes  Yes 

Leave  0  for  general  use 

IGLMAX  Number  of  sweeps. 

Low  Reynolds  number 
flows  need  more  sweeps 
to  converge . 


25 


IPRNT 


Print  flag  to  control 
output  print. 


1 


Summary  print  for  the 
last  two  sweeps 


=  0  Summary  print  for  each 
sweep . 

=  2  Whole  print  for  all 
sweeps . 


[  Control  Line  2  I  FORMAT 


Column 

Name 

1  -  10 

RL 

11-20 

XTRL 

21  -30 

XTRU 

31-40 

ALPO 

I  Control  Line  3  I 

Column  Name 

1  -5  MPTS 


I  Coordinate  Data  Line  i 


(  4  E10 . 0  ) 

Explanat ion 

Chord  Reynolds  number 

Fixed  transition  location 
for  lower  surface 

Fixed  transition  location 
for  upper  surface  when  the 
stagnation  point  is  above 
the  leading  edge,  XTRL  may 
be  positive.  When  the 
stagnation  point  is  beneath 
the  leading  edge,  XTRU  may 
be  negative. 

Angle  of  attack,  in  degree 


FORMAT  (15) 

Explanation 

Number  of  airfoil 
coord ina  tes . 


FORMAT  ( 2F  10.0) 


Input  dimensionless  airfoil  coordinates  as  two  columns  in 
one  row  (  x/c,  y/c).  The  order  is  trailing  edge  — *  lower 

surface  — »  leading  edge  - *  upper  surface  - *  trailing 

edge.  Thus,  trailing  edge  is  input  twice. 


B.  OUTPUT  DESCRIPTION 

1 .  Introduction 

To  take  the  appropriate  results  from  the  output, 
user  must  check  the  convergence.  This  can  be  done  by 
comparing  the  convergence  indicators,  lift  coefficient  and 
displacement  thickness  at  trailing  edge,  when  the 
computation  is  completed  successfully  for  the  given  sweeps. 
If  the  convergence  indicators  show  steady  values  over  some 
sweeps  i.e.  each  differ  nee  of  the  convergence  indicators  is 
less  than  1%,  the  results  are  considered  as  converged. 

Sometimes,  the  user  may  experience  failure  in 
computations  at  low  Reynolds  numbers  and  high  angles  of 
attack  due  to  numerical  breakdowns.  These  breakdowns  take 
place  when  Newton  iteration  does  not  converge  near  to  the 
trailing  edge  on  the  upper  surface  or  the  computation  is 
terminated  by  Fortran  error  with  this  message: 

IFY  2511  SSQRT;  ARG  =  argument,  LESS  THAN  ZERO 

Also,  extensive  flow  separation  cause  numerical 
breakdowns  of  the  boundary  layer  computation.  These 
unrealistically  large  regions  of  separation  at  low  Reynolds 
numbers  can  be  reduced  by  changing  the  transitional  flow 
mode  1  . 

2 .  Detailed  Output  Data  Description 


The  output  of  the  computer  program  for  IPRNT 


and  0  can  be  divided  into  three  parts  as  follows: 

1.  Input  data  and  inviscid  lift  coefficient 

2.  Alternating  boundary  layer  parameters  for  each 
surface  and  wake. 

3.  Inviscid  and  viscous  pressure  distributions 

The  following  is  the  list  of  computed  boundary  layer 
parameters  in  the  order  in  which  they  appear  in  the  output 


print 


Name 

in  print 


Name 

in  program 


Defini t ion 


Explanation 


X  (  NX  ) 

X(  I  ) 

J_ 

Dimensionless  sur¬ 

c 

face  distance  from 
the  stagnation  point 

X/C( NX ) 

XC  (  I  ) 

X 

c 

Dimensionless  chord- 
wise  distance  from 
leading  edge 

V( 1,NX) 

VWL ( I ) 

/  d  "  f 

\  Dimensionless  wall 

\  d  V 

/w  shear  stress  para¬ 

meter. 

CF 

CF  (  I  ) 

2 

Local  skin  friction 

coefficient 

DELST 

DLS ( I ) 

J  * 

Dimensionless  dis¬ 

c 

placement  thickness. 

l(‘ 

u  )  dy 

U«  ' 

THETA 

THETA ( I ) 

_0 

Dimensionless 

C 

,00 

momentum  thickness 

e  = 

u  f. 

u  \  dy 

J  Ue  i 

Q 

U*  / 

UE(NX) 

UE(  I  ) 

Ue  I 

Dimensionless  exter¬ 

u. 

nal  velocity 

computed  from  the 
inviscid  method  with 
viscous  effect. 
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W ( NP, NX )  WNP(I) 


0 ( NX )  D ( KX ) 

DB(NX)  DB(KX) 


UeV 

U. 


Ue  v  £*Jrl 

U  co  c 


Dimensionless  exter¬ 
nal  velocity 
computed  from  the 
interactive  boundary 
layer  method. 

Product  of  velocity 
and  displacement 
thickness  in  dimen¬ 
sionless  for m  f o r 
the  current  sweep 
( D ) ,  from  the  pre¬ 
vious  sweep  ( DB  )  . 

D  and  DB  are  printed 
in  the  output  as  up¬ 
dated  values  through 
the  relaxation 
procedure 


IT 


(N*)  r  £,j  (**)[  I  d-  ( 


W  (np,  \x) 
U 


ITN(I)  Iteration  count  for 

the  convergence  of 
solutions  at  a  given 
NX-station . 


NP  NNP(I) 


Number  of  points 
across  the  boundary 
layer  . 


C.  HOW  TO  CHANGE  THE  ORIGINAL  PROGRAM 

Cebeci's  computer  program  is  written  entirely  in  FORTRAN 
and  consists  of  the  following  eight  FORTRAN  files: 


FILE 

1 

FORTRAN 

A 1 

FILE 

2 

FORTRAN 

A  1 

FILE 

3 

FORTRAN 

A 1 

FILE 

4 

FORTRAN 

A  1 

FORTftfrtt- 

FILE  6 


FORTRAN  A 1 


FILE  7  FORTRAN  A1 


FILE 

8 

FORTRAN 

A 1 

FILE 

9 

FORTRAN 

A 1 

Also,  the  user  needs  the  following  six  JCL  files  to  run 
changed  program: 


QUEST  1 

JCL 

A 1 

ALLPDS 

JCL 

A 1 

QUEST  3 

JCL 

A 1 

STOPDS 

JCL 

A 1 

LOADMOD 

JCL 

A 1 

INTAPCLC 

JCL 

A 1 

"QUEST  1"  lists  all  data  set  on  MSS  (Mass  Storage 

System)  that  belong  to  student  user  ID 
number  and  in  group  PUB4B  (public  disk, 
volume),  which  deletes  the  data  set  180 
days  after  creation. 

"ALLPDS"  allocates  space  fo^  a  PDS  (partitioned 

Data  Set)  on  PUB4B  prior  to  reading  a 
data  set  into  it,  using  the  utility 
program  IEFBR  14  which  pre-al locates  or 
deletes  data  sets. 

"QUEST  3"  lists  the  members  of  a  PDS  and  calculate 

the  remaining  available  space  within  the 
data  set. 

"S70PDS"  places  a  FORTRAN  source  file  in  a  PDS. 


The  SYSUT2  DD  statement  describes  the  out¬ 
put  data  set  and  the  SYSUT1  DD  statement 
describes  the  input  data  set. 

"LOADMOD"  creates  a  load  module. 

"  I NTAPCLG "  compiles  scarce,  cede,  loads  qnd  O-Ccof^s  Text  file 

A  list  of  all  JCL  files  is  given  at  the  end  of  this  chapter 
and  more  information  is  in  the  User's  Guide  to  MVS  at  NPS 
(  1986  )  . 

The  following  is  the  procedure  to  change  the  program 
and  run  it: 

Step  1 

1.  Make  changes  in  FORTRAN  files  as  required.  Remember 
those  eight  FORTRAN  files  compose  one  program. 

Thus,  the  user  must  check  all  files  which  are 
related  to  changes. 

2.  Review  and  update  FORTRAN  errors,  if  necessary. 

Step  2 

1.  Change  job  name  "xxxxxx",  user  ID  number  "nnnn"  and 
library  name  "yyy"  with  user's  own  in  all  JCL  files 
to  make  them  be  the  user's. 

2.  Submit  four  JCL  files: 

"QUEST1 "  " ALLPDS"  "QUEST  3"  "LOADMOD". 

Step  3 

Do  the  following  for  all  FORTRAN  files  to  submit 
them  to  user's  library. 

1.  Open  " STOPDS"  JCL  file. 

2.  Change  the  FORTRAN  file  number  indicated  by  n  on  the 
sixth  line. 


/ /SYSUT2  DD  DISP 


YYYLIB  (INTAIR  n) 


3.  Insert  the  FORTRAN  file  between 

"//SYSUT1  DD,  *"  and 

4.  Submit  ''STOPDS'’  JCL  file. 

5.  Check  the  result  by  any  error  message  at  the  very 
top  and  bottom  part  and  condition  code. 

6.  Update  errors,  if  necessary. 

Step  4 

Do  the  following  to  run  the  changed  program. 

1.  Open  "INTAPCLG"  JCL  file. 

2.  Insert  the  input  data  between 

"//GO.  SYSIN  DD  *  "  and  "/*". 

3.  Submit  "INTAPCLG"  JCL  file 

From  the  second  change,  do  steps  1,  3  and  4. 

D.  APPLICATION  OF  CEBECI’S  PROGRAM 

Cebeci's  interactive  program  was  applied  to  a  single 
airfoil,  FX  63  -137,  at  three  Reynolds  numbers.  The  print 
flag  is  0  and  the  number  of  airfoil  coordinates  is  49  for 
ail  cases.  The  results  are  compared  with  experimental  data 
and  the  turbulence  model  has  been  changed  to  get  better 
results.  The  experimental  results  are  taken  from  Reference 
9  . 

Figure  4.2  and  4.3  show  the  comparison  of  lift  and  drag 
coefficients  for  three  Reynolds  numbers.  Drag  coefficients 
are  computed  at  the  trailing  edge.  It  is  seen  that  the 
results  are  closer  to  the  experimental  data  as  the  Reynolds 


number  increases. 


Figures  4.4,  4.5  and  4.6  show  the  skin  friction 
coefficients  on  the  upper  and  lower  surfaces  as  a  function 
of  angle  of  attack  for  the  same  three  Reynolds  numbers. 
Legends  indicate  the  angle  of  attack.  As  the  Reynolds 
number  increases,  the  length  of  the  separation  bubble 
decreases  on  both  surfaces  for  the  same  angle  of  attack. 
Especially,  Figure  4.4,  which  has  a  different  label  of 
y-axis  from  the  other  figures,  indicates  an  unrealistic  flow 
on  the  upper  surface  at  low  Reynolds  number.  In  Figure  4.5 
and  4.6,  however,  the  flow  on  the  upper  surface  is  separated 
from  about  30%  chord  at  ail  angles  of  attack. 

Figures  4.7,  4.8  and  4.9  present  how  the  displacement 
thickness  varies  according  to  the  angle  of  attack  and  the 
Reynolds  number.  As  the  Reynolds  number  increases,  the 
displacement  thickness  decreases  on  both  surfaces.  As  the 
angle  of  attack  increases,  it  also  increases  on  the  upper 
surface,  but  decreases  on  the  lower  surface. 

Table  4.1  represents  the  computed  and  fixed  transition 
locations  which  are  used  in  Figure'  4.10  4.13.  The 

computed  transition  locations  XTRL  =  C. 288922  at  ALPD  =  8 
degrees  and  XTRL  =  0.206866  at  ALPD  =  10  degrees  were 
printed  incorrectly.  In  this  case,  the  correct  transition 
locations  can  be  obtained  by  the  following  procedure: 

1.  Calculate  the  difference  of  XTRL's  at  ALPD  =  4 
degrees  and  6  degrees. 

2.  Add  the  difference  to  XTRL  of  ALPD  -  6  degrees  and 


take  the  result  as  a  temporary  XTRL  for  ALPD  = 

8  degrees. 

Fix  transition  location  with  the  temporary  XTRL 


3  . 

4.  Run  the  program  to  print  XC(I),  GAMTR(I)  from 
SUBROUTINE  OUTPUT  in  FILE  3  FORTRAN  Al. 

5.  Take  XC ( I ) ,  where  the  value  of  GAMTR(I)  is  finally 
equal  to  zero,  as  the  correct  transition  location. 

6.  From  the  same  procedures  with  XTRL's  at  ALPD  =  6 
degrees  and  8  degrees,  obtain  the  correct  transition 
location  for  ALPD  =  10  degrees. 

Figures  4.10-^4.13  show  the  effect  of  variations  in  the 

empirical  constant,  Gy*  ,  in  turbulence  model  at  Re  -  0.28 

6 

x  10  In  the  legend  of  Figures  4. 1C  and  4.11,  the  three 

numbers  are  the  values  of  the  empirical  constant,  C  means 
that  the  transition  location  is  computed  by  Eq .  (4.1)  and  F 

means  that  it  is  fixed.  The  lift  curve  is  closer  to  the 
experimental  data  as  the  empirical  constant  decreases,  but 
the  drag  curve  is  closest  when  the  empirical  constant  is 
120  . 

The  reason  why  we  can  obtain  better  results  by  reducing 
the  empirical  constant  is  shown  in  Figure  4.12  and  4.13. 
GAMMATR  (  ftr)  is  the  inter  mi ttency  factor  which  is  a 
function  of  the  x-coordinate  with  values  0.0  at  the 
beginning  point  of  transition  and  1.0  at  the  ending  point  of 
transition.  The  inter m i t tency  factor  makes  it  possible  to 
avoid  a  sudden  transition  from  laminar  to  turbulent  by 
smoothing  out  the  step-shaped  change  of  viscosity  from 
kinematic  to  eddy.  The  relation  between  the  empirical 


_  j'^vy.V.V Jjf 


constant  and  the  inter mi t tency  factor  is  given  as: 

r  j  -  i .  3*  ,x 

/<r  =  1  -  exp  -  (x  -  xtr)  Ue  P.^x  1  dx 

L  Grtr  ^  JXtr 

where,  x^^.  denotes  the  beginning  point  of  transition. 

The  transition  length  is  mainly  determined  by  the 

empirical  constant  Gy  which  is  set  as  1203  in  CeLeci's 

ytr 

original  program.  Decreasing  the  value  of  the  empirical 
constant  reduces  the  transition  length.  In  this  way,  the 
unrealistic  flow  shown  in  Figure  4.4  can  be  avoided  and  more 
reasonable  results,  which  are  closer  to  the  experimental 
data,  can  be  obtained.  In  Figure  4.12  and  4.13,  the  legend 
indicates  the  angle  of  attack,  and  each  of  four  line 

patterns  is  used  twice  for  two  angles  of  attack  (e.g.  -  ; 

4  degrees  and  4  degrees).  The  direction  of  curves,  as  the 
angle  of  attack  increases,  is  from  right  to  left  at  the 
upper  surface  ^Figure  4.12)  and  from  left  to  right  at  the 


lower  s ur face 


(  F iyur e  4.13). 


List  of  JCL  Files 


xxxxxxx  ALLPDS  JCL  A1  xxxxxxx 


//XXXXXX  JOB  (NNNN. 9999), 'ALLOCATE  PDS’,CLA5S-A 
//  XMA I N  0RG=NPGVM1 .NNNNP 
//  EXEC  PGM  -  I EFBR 1 4 

//SYSPRINT  DD  SYSOUT^A 

// DD1  DD  UNI T  =33J0V,M3VGP-PUB4B. DISP=(NEH,CATLG, DELETE) , 
//  SPACE- (CYL ,(4,4,61), DSN=MSS . SNNNN . YYYLIB 

// 


xxxxxxx  QUEST  1  JCL  A1  xxxxxxx 


//XXXXXX  JOB  (NNNN, 9499), 'QUESTION* ,  C  L  A  S  S  -  A 

//xMAIN  0RG=NPGVM1 .NNNNP 

//  EXEC  PGM= I DCAMS 

//SYSPRINT  DD  SYS0UT=A 

//SYSIN  DD  x 

LISTDSET  GROUPC  PUB4B )  L EVEL ( MSS . SNNNN ) 

/  x 

// 


xxxxxxx  QVJE3T3  JCL  A1  xxxxxxx 


//XXXXXX  JOB  (NNNN, 9999), 'LIST' ,CLASS=A 
//XMAIN  0RG=NPGVM1 . NNNNP 
//  EXEC  PGM^IEHLIST 
//SYSPRINT  DD  5YS0UT=A 

// DD1  DD  UNIT=3330V, VOL =SER=M50005 , DISP=5HR 
//SYSIN  DD  x 

L ISTVTOC  FORMAT, VOL=3330V-M50005, DSN AME:M5S . SNNNN. YYYLIB 
L ISTPDS  V0L-3330V=MS0005,DSNAME=MSS . SNNNN. YYYLIB 
/  x 

// 


xxxxxxx  STOPDS  JCL  A1  xxxxxxx 


//XXXXXX  JOB  (NNNN, 9999), 'PLACE  PDS' 

//•MAIN  0RG-NPGVM1 .NNNNP 
//  EXEC  PGM=  I EBGENER 

//SYSPRINT  DD  SYSOUT^A 
//SYSIN  DD  DUMMY 

//SYSUT2  DD  DISP  =  ( OLD, KEEP), DSNAME;MSS . SNNNN. YYYL I B ( I  NT  A I R_ ) , 
//  DCB-(RECFM-FB, LRECL  -  80, Bi.KSIZE-8000  ) 

//SYSUT1  DD  * 

/  x 

// 


hQ 


'•-.1 


.  .  .  .  .  ,  .  O  V*  v'  \  1 


xxxxxxx 


//XXXXX  JOB  CNNNN, 9999), 'CREATE  L OADMODUL E ' , CL ASS=C 
//XMAIN  ORG-NPGVM1 . NNNNP 

//  EXEC  FORTVCL.PARM. FORT= 'LVL(77),NOS,MOX, NOMAP ' 

//FORT . SYSPRINT  DD  DUMMY 

//FORT . SYS  I N  PD  DI SP  =  SHR , DSN  =  MSS . SHNNN . YYYL I B ( I  NT A I R1 ) 

//  DD  DI SP=S HR , DSN = MSS . SNNNN . YYYL IBCINTAIR2) 

//  DD  DI SP=5HR , DSN-MSS . SNNNN .  YYYL  IB(  INTAIR3) 

//  DD  DI3P=SHR,  DSM-M3S  .  SNMNtJ .  YYYL  I  B  C  INTAIR4) 

//  DD  DI 5P=SHR , DSN-MSS . SNNNN .  YYYL  IBCINTAIR6) 

//  DD  D I SP  =  5HR ,  DSN -MSS , SNNNN . YYYL IBC INTAIR7  ) 

//  DD  DISP=SHR,  DSN=MSS  .  SNNNN  .  YYYL  IB(  INTAIR3) 

//  DD  DISP=3HR, DS Nr MSS . SNNNN . YYYL IBCINTAIR9) 

//LKED.  SYSLMOD  DD  DI SP -SHR , DSN-MSS . SNNNN .  YYYL  IBC INTCAS) 
//LKED.SYSIN  DD  DSM-NULLFILE 
// 


x  x  x  x  x  x  x  INTAPCLG  JCL  A1  xxxxxxx 


//XXXXXX  JOB  (NNNN, 9999), 'COMP  LOAD  AND  GO',CLASS=G 
//*MAIN  ORG=NPGVMl . HNNNP 

//  EXEC  FORTVCLG, PARM . FORT= • LVL < 77 ) , NOS, NOX, NOMAP' , 

//  REGION. GO=1280K 

//FORT . SYSPRINT  DD  DUMMY 

//FORT . SYGIN  DD  DI SP =SHR , DSN=MSS . SNNNN . YYYL I B ( INTAIRI ) 

//  DD  DISPrSHR, DSN =MSS. SNNNN. YYYL IBC I NTAIR2) 

//  DD  DISP=SHR, DSN=MSS. SNNNN. YYYLIBC INTAIR5) 

//  DD  DISP=SHR, DSN=MSS. SNNNN. YYYLIBC INTAIR4) 

//  DD  DISP  =  SHR,DSN  =  MSS. SNNNN . YYYL IBC INTAIR6  ) 

//  DD  DI SP=SHR , DSN : MSS . SNNNN . YYYL IBC INTAIR7 ) 

//  DD  DISP=SHR, DSN -MSS. SNNNN. YYYLIBC I NTAIR8) 

//  DD  DISP=SHR,DSN=MSS. SNNNN . YYYLIBC INTAIR9 ) 

//LKED. SYSPRINT  DD  DUMMY 

//LKED.SYSIN  DD  D5N-NULLFILE 

//GO . FTOlFOOl  DD  UNI T  =  SYSDA , SPACE-CCYL , C 1 , 1 ) ) 

//GO  .  FT02F001  DD  UNI T  =  SYSDA , SPACE  =  C CYL , C 1 , 1 )  ) 

//GO . FT03F001  DD  UNIT=SYSDA,SPACEICCYL, C 1 , I ) ) 

//GO . FT04F001  DD  UNIT  =  SYSDA, SPACE  =  CCYL , ( 1 , I  )  ) 

//GO . FT99F001  DD  UNIT  =  SY5DA , SPACEZCCYL , C I , 1 ) ) 
//GO.FT08F001  DD  UNI T =SYSDA , SPACE = C CYL , C 1 , 1 ) ) 

//GO . FT09F001  DD  UNI T=SYSDA , SPACE= C CYL , ( I , 1 ) ) 

//GO . FTI0F001  DD  UNIT  -  SYSDA, SPACE=CCYL , C 1 , 1 )  ) 
//GO.FT11FOOI  DD  UN  I T =3Y5 DA , SPACE  =  C CYL , C 1 , 1 ) ) 

//GO . FT06F00I  DD  SY50UT=A 
//GO.SYSIN  DD  x 
/  x 
// 
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Cebeci's  interactive  computer  program  was  applied  to 
the  Wor tmann-Al thaus  FX  63-137  airfoil  to  show  the 
capability  of  strong  v iscous / i nv  i  sc  id  interaction  methods  to 
predict  airfoil  flows  at  low  Reynolds  numbers. 

From  the  comparisons  with  the  experimental  data,  it 
was  confirmed  that  the  results  are  closer  to  the 
experimental  data  as  the  Reynolds  number  increases. 

Also,  much  better  results  were  obtained  by  decreasing 

the  empirical  constant  G^- 

Mr 

Therefore,  it  was  concluded  that  the  boundary  layer 
transition  model  has  an  important  influence  on  the 
predictive  capability  of  viscous/ inviscid  interaction 
methods . 
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